function fI = SpectralInterp(varargin)

%Interpolates from spectral grids (either fourier or cheb) to arbitrary
%grid.  Fourier grid is assumed to be x = -pi + 2*pi*(0:N-1)/N and cheb
%grid is assumed to be x = cos(pi*(2*(1:N)-1)/(2*N)).  Input format is
%(f,xI,kind,dim).  dim input is set to 1 if dim is not inputed. kind is a
%string which is either 'fourier' or 'cheb'


f    = varargin{1};
xI   = varargin{2}(:);
kind = varargin{3};
if length(varargin) == 4
    dim = varargin{4};
else
    dim = 1;
end

sizef = size(f);
N     = sizef(dim);
NI    = length(xI);
Ndim  = ndims(f);

tol = 1e-2;

if strcmp(kind,'cheb')
    
    m  = (1:N)'; m = m(:,ones(1,NI));
    x  = cos(pi*(2*m-1)/(2*N));
    xI = xI(:,ones(1,N)); xI = xI';  
    I  = find(abs(x-xI)<tol); 
     
    Cmn = (-1).^m.*cos(N*acos(xI)).*sin((pi-2*pi*m)/(2*N))./(N*(xI-cos((pi-2*pi*m)/(2*N)))); %this will be badly behaved near cheby points
    
    if numel(I) > 0
        Cmn(I) = 0;
        
        for k = 0:N-1
            if k == 0
                c = 1;
            else 
                c = 2;
            end
            Cmn(I) = Cmn(I) + (c/N*cos(k*(pi*(2*m(I)-1)/(2*N)))).*cos(k*acos(xI(I)));
        end
    end
    
    F = multiprod(Cmn,f,1,dim);
    
    order = [2:dim,1,dim+1:Ndim];
    
    fI = permute(F,order);

    
elseif strcmp(kind,'cheb2')
    N = N - 1;
    m  = (0:N)'; m = m(:,ones(1,NI));
    x = cos(pi*m/N);
    xI = xI(:,ones(1,N+1)); xI = xI';  
    I  = find(abs(x-xI)<tol); 
    
    cm = ones(size(m)); cm(abs(m) == N) = 2;  cm(abs(m) == 0) = 2;
     
    Cmn = -(-1).^(m).*sqrt(1-xI.^2).*sin(N*acos(xI))./(xI - cos(pi*m/N))./cm./N; %this will be badly behaved near cheby points
    
    if numel(I) > 0
        Cmn(I) = 0;
        
        for k = 0:N
            if k == 0 || k == N
                p = 2;
            else 
                p = 1;
            end
            Cmn(I) = Cmn(I) + cos(k*pi*m(I)/N).*cos(k*acos(xI(I)))./cm(I)/p*2/N;
        end
    end    
            
    F = multiprod(Cmn,f,1,dim);
    
    order = [2:dim,1,dim+1:Ndim];
    
    fI = permute(F,order);

    
elseif strcmp(kind,'fourier')
    
    x = -pi + 2*pi*(0:N-1)'/N; x = x(:,ones(1,NI));
    xI = xI(:,ones(1,N)); xI = xI'; 
    I  = find(abs(x-xI)<tol);
    
    Cmn = 1/N*csc((xI - x)/2).*sin(N/2*(xI - x)); %this will be badly behaved near fourier points
    
    if numel(I) > 0
        M = (N-1)/2;
        Cmn(I) = 0;
        
        for k = 0:M
            if k == 0
                c = 1;
            else 
                c = 2;
            end
            Cmn(I) = Cmn(I) + c/N*cos(k*(xI(I) - x(I)));
        end
    end

    F = multiprod(Cmn,f,1,dim);
    
    order = [2:dim,1,dim+1:Ndim];
        
    fI = permute(F,order);

else
    error('unknown kind')
end

function C = multiprod(A,B,dimA,dimB)

%multidimensional array product
%C(a1,a2,...aN,b1,b2...bM) = A(a1,a2,...k,...aN)*B(b1,b2...k,...bM)

sizeA = size(A);
sizeB = size(B);
NdimA = sizeA(dimA); %number of elements in dimA
NdimB = sizeB(dimB); %number of elements in dimB

if NdimA ~= NdimB
    error('contraction dimensions must agree')
end

%permute dimA of A to last dim and dimB to first dim
orderA = 1:length(sizeA); orderA(dimA) = []; orderA = [orderA, dimA];
orderB = 1:length(sizeB); orderB(dimB) = []; orderB = [dimB,orderB];
AA = permute(A,orderA);
BB = permute(B,orderB);

%reshape into matricies
NA = numel(A)/NdimA;
NB = numel(B)/NdimB;
A2 = reshape(AA,NA,NdimA);
B2 = reshape(BB,NdimB,NB);

s1 = sizeA; s1(dimA) = [];
s2 = sizeB; s2(dimB) = [];
sizeC = [s1,s2];

C = reshape(A2*B2,sizeC);
